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Abstract 

Exploratory variational pseudopotential density functional calculations are performed for the 
electronic properties of many-electron systems in the 3D cartesian coordinate grid (CCG). The 
atom-centered localized gaussian basis set, electronic density and the two-body potentials are set 
up in the 3D cubic box. The classical Hartree potential is calculated accurately and efficiently 
through a Fourier convolution technique. As a first step, simple local density functionals of homo- 
geneous electron gas are used for the exchange-correlation potential, while Hay-Wadt-type effective 
core potentials are employed to eliminate the core electrons. No auxiliary basis set is invoked. 
Preliminary illustrative calculations on total energies, individual energy components, eigenvalues, 
potential energy curves, ionization energies, atomization energies of a set of 12 molecules show 
excellent agreement with the corresponding reference values of atom-centered grid as well as the 
grid-free calculation. Results for 3 atoms are also given. Combination of CCG and the convolution 
procedure used for classical Coulomb potential can provide reasonably accurate and reliable results 
for many-electron systems. 
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I. INTRODUCTION 



In the past few decades, density functional theory (DFT) [l|, |2| based calculations have 
played a pivotal role in our understanding of the electronic structure, properties and dy- 
namics of many-electron systems such as atoms, molecules and solids. This overwhelming 
success is chiefly due to their ability to offer accurate results and to account for the electron 
correlation effects in a transparent and tractable manner. Accordingly, a large number of 
elegant and efficient approaches are available today for the treatment of these systems cov- 
ering a broad range of approximations, accuracy, computational and algorithmic schemes. 
Consequently, this has been a very fertile area of research in the recent years and continues 
to remain at the forefront of modern research. 

Leaving aside a few attempts such as the high-order finite difference method j^, where 
without the explicit use of basis set, the discretized Kohn-Sham (KS) equation is directly 
solved on real-space grid, a vast majority of today's all-electron or pseudopotential DFT 
methodologies instead employ some suitable basis set. Some of the notable ones include 
the plane waves (PW) or atom-centered localized basis sets such as Slater type orbitals 
(STO), Gaussian type orbitals (GTO), numerical radial functions, linear muffin-tin orbitals, 
etc. Among these, the GTO basis sets have gained more popularity over the others for 
nonperiodic systems such as molecules and clusters due to the convenient analytic routes 
they provide for the relevant matrix elements involving multiple centers. On a similar 
ground, PW basis sets are the preferred options for periodic systems. A gaussian and 
augmented-plane-wave DFT approach has also been reported g, |7| where the KS molecular 
orbitals (MO) and the electronic charge densities are expanded in the gaussian basis and an 
augmented PW basis sets respectively. 

Nowadays, linear combination of GTO based DFT calculations have become an invaluable 
and routine tool for quantum chemists, where the KS MOs are typically expanded in terms 
of the contracted gaussian functions centered on the atoms^, ISj. Furthermore, the electron 
density as well as the exchange-correlation (XC) potentials are also similarly expanded 
in terms of a finite number of auxiliary gaussian type basis functions. This obviates the 
necessity to compute the expensive four-centered integrals, making it an A^^ process. 

In the grid-based methods, typically the multi-center molecular integrals are decomposed 
into monocentric atomic subintegrals using some weighting scheme and then each of these 
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a. 



latter are computed numerically on the respective atomic grid [8[. The radial integrals 
a.e performed though a vanety of technique. Gau.-Cheby=hev ,uad.atu.e of second 

kind [8||, Chebyshev quadratures of first and second kind in conjunction with several map- 
ping schemes bll, simple Gauss quadrature 10|, transformation based on Euler-MacLaurin 

n n n 

formula [11|, 112|, numerical quadratures [13|, etc. Angular integrations often use the octa- 



hedral grid developed by Lebedev [14] or the Lobatto approaches [9[. Automatic numerical 
integrator generating an adaptive molecular grid depending on size and shell structure of the 
given basis set, has also been developed 15|. Recently, there has been interest in the real- 
space cartesian grid as well. For example in the Fourier Transform Coulomb method jl^, 
molecular integrations are performed by dividing the gaussian shell pairs into "smooth" and 
"sharp" categories based on their exponents. Of late, a multiresolution technique, combining 
the atom-centered and cubic cartesian grid, with divided difference interpolation pla yin g the 
role of a communicator between the two, have been shown to be quite efficient 17|, Il8 |. 

In this article, we make a detailed systematic analysis of the performance and relevance 
of the cartesian coordinate grids(CCG) in the context of molecular DFT calculations. As 
our results in the following sections show, this is indeed capable of producing fairly accurate 
and physically meaningful results, at least for small molecules. We use the usual GTO- 
type linear expansion of the KS MOs; however no additional auxiliary basis set is utilized 
to express the charge density or the XC potential. The basis functions, the MOs, the 
electron density as well as the various two-electron potentials are directly set up in the 3D 
real CCG. The classical Hartree potential is obtained accurately and efficiently through a 



Fourier convolution method involving a set of FFT-inverse FFT pair 2l|, |2^. Analytical 



one-electron ab initio effective core potentials, consisting of a sum of gaussian functions 
are employed to represent the core electrons while energy-optimized truncated gaussian 
bases are chosen for the valence electrons. Here we restrict ourselves to the local density 
functionals of homogeneous electron gas to incorporate the XC effects, while more accurate, 
sophisticated functionals may be considered in future. The KS matrix eigenvalue equation is 
solved in the usual self consistent manner to obtain the KS orbitals and eigenvalues. Results 
are given for 12 representative molecules (5 diatomics and 7 polyatomics) and 3 atoms. 
In order to assess the accuracy and reliability, first we make a systematic investigation on 
the convergence of total energies, individual energy components and the integrated electron 
density, in various grid representation through a comparison with the reference theoretical 
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results. Then we compare the eigenvalues, atomization energies and the ionization energies 
with the hterature data. Potential energy curves are given for two diatomics (HCl and CI2). 
The article in organized as follows: Section II gives the methods of computation. Section 
III presents a discussion on the results obtained while a few concluding remarks are made 
on the future and prospect in Section IV. 



II. METHOD OF CALCULATION 



We are interested in the single-particle KS equation for the ground state of a many- 
electron system under the influence of pseudopotential, which can be written as (atomic 
units employed unless otherwise stated), 

1. 



-^V^ + <on(r) + VH[p\{r) + vxc[p\{r) 



V'i(r) = eiijjiir) (1) 
where vf^^ denotes the ionic pseudopotential for the system. 



<on(r) = E^ka(r-Ra) (2) 

with f fo„ „ signifying the ion-core pseudopotential associated with atom A, situated at Ra. 
The Hartree potential t'H[p](r) describes the classical electrostatic interactions among va- 
lence electrons while the XC potential vxc[p\{'>^) represents the non-classical part of the 
Hamiltonian. {■0f(r)},(j — a or corresponds to a set of N occupied orthonormal MOs 
and the electronic density is given by the expression, 

p(r) =p«(r)+pV) = E/riC(r)r + E/fl^f(r)r (3) 

i i 

where the /j's denote the occupation numbers. 

The basis functions and the MOs are directly built on the real uniform 3D cartesian grid 
simulating a cubic box, 

n^ri^ {i-l)hr, i ^ 1,2,- ■■ ,Nr; for re{x,y,z} (4) 

where hr and Nr denote the grid spacing and number of points in the grid respectively. The 
electron density in the grid is given by, 

Pi^g) = E Pi^uXi,{rg)xArg) (5) 

HI/ 



where P^^p signifies an element of tlie density matrix and tlie set {x^(r)} contains tlie con- 
tracted gaussian functions centered on tlie atoms. Any quantity including those in the above 
will be represented in the discretized grid with a subscript "g" . 

The standard matrix form of the KS equation reads (for a system having K basis func- 
tions), 

F^^C = SCe (6) 

where C is a X x X square matrix containing the expansion coefficients C^j and e is the 
diagonal matrix of the orbital energies ej. S corresponds to the overlap matrix and F^'^ 
denotes the total KS matrix including the core Hamiltonian and the effective KS potential, 

^/f/ = / X,ij) [/i— + VHxcij)] x.(r) dv = + {x,{r)\vHxc{r)\xAr)) (7) 

The first term in the right hand side accounts for the one-electron energies in the Hamil- 
tonian. The overlap, kinetic and nuclear-electron attraction integrals are identical to those 
obtained in the gaussian basis-set based HF methods and we have used the standard re- 
cursion algorithms 



231, 



2J] for their evaluation. Significant progress has been made in 



the development of rigorous ab initio effective core potentia 
ploy the angular-momentum dependent form as proposed by 



s and in this work we em- 
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261]. While construction 



of the pseudopotential matrix elements in gaussian basis sets is currently in progress, in 
the present implementation, they are imported from the widely used quantum chemistry 
program GAMESS output (with printing option S)^. 

For finite systems, the simplest and perhaps the most crude way to calculate vh is hj 
direct numerical integration, which is feasible for only small systems. However the most 
widely used approach is to solve the corresponding Poisson equation. Recently however, in 
the literature, the conventional Fourier convolution method and some of its variants have 
been shown to be quite accurate and efficient in the context of molecular modeling. In 



particular, here we adopt the approach as put forth in [21 



22|. 



p(k,) = FFT{p(r,)} 
VH{r,) = FFT-i{^;|,(k,) p(k,)} (8) 

Here p(kg) and vj^ikg) denote the Fourier integrals of density and the Coulomb interaction 
kernel in the grid. It may be noted that while /9(kg) can be easily obtained from the discrete 
Fourier transform of its real-space values by standard FFT, the calculation of the latter is 



nontrivial and requires caution due to the presence of singularity in the real space. This is 
tackled by decomposing the kernel into long- and short-range terms: 

. , , erf(ar) erfc(ar) r / n r / \ 
^nir,) = + ^ Vk^^Jr,) + Vk^^^Jr,) (9) 

where erf(x) and erfc(x) denote the error function and its complements respectively. The 
Fourier integral of the short-range part is calculated analytically whereas the long-range 
interaction is directly obtained from the FFT of the real-space values. It may be mentioned 
that while the above FFT-based Ewald type summation method scales as In A^, two other 
Poisson-solvers have gained popularity in the context of large-scale electronic structure cal- 
culation within the KS DFT framework, which are quite efficient and scale linearly. In 
the fast multipole method (FMM) type approaches, the near-field contributions are tackled 
explicitly, whereas the far-field is treated through a clustering of the spatial cells and rep- 
resenting the field with a multipole expansion. Another route employs the highly efficient 
multigrid method within the real-space formalism such as a finite-difference or finite-element 
scheme. A review of the various available techniques could be found in 



One of the most critical issues in any DFT calculation is the choice of appropriate XC 
functionals, the exact form of which is unknown as yet and must be approximated. In 
the literature, an enormous number of functionals with varying complexity, property and 
accuracy have been published; the present work employs the simple local XC functionals of a 



homogeneous electron gas (formula V of ref. 291]). In the absence of any analytical method. 



the corresponding two-electron matrix elements can be either calculated numerica. 
grid or be fitted by an auxiliary set of gaussian functions as suggested by 
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ly on the 



employed in some of the existing DFT codes 



3l| and 



This work employs the direct numerical 



integration on the CCG to obtain these matrix elements, viz., 

(X/.(r<;) \vHXc{rg)\xu{rg)) = Khyh, ^ X/.(r<;) ^Hxcirg) Xu{rg) (10) 

9 

The algebraic eigenvalue equation is solved in the usual self-consistent manner and the total 
energy of the system can be obtained as a sum of the various components in the standard 
way. Three convergence criteria were imposed in this work viz., (i) the electronic energy 
differences between two successive iterations is below certain threshold (ii) the maximum 
absolute deviation in the potential is less than the specified tolerance limit and (iii) the 
standard deviation in a density matrix element remains within a prescribed threshold. For 



TABLE I: Variation of all the energy components and total number of electrons with respect to 
the grid parameters for CI2 along with the reference values at R = 4.20 a.u. All quantities are in 



a.u. 





Nr = 


32 




= 64 






128 


Nr = 256 


Ref. [27] 


Set 


A 

A 


B 


C 


D 


E 


: 


G 


H 




hr 


0.3 


0.4 


0.2 


0.3 


0.4 


0.1 


0.2 


0.1 




(T) 


11.00750 


11.17919 


11.18733 


11.07195 


11.06448 


"1 "1 "1 "1 

11.18701 


11.07244 


11.07244 


11.07320 


\ t 1 


-83.43381 


-83.68501 


-83.70054 


-83.45722 


-83.44290 


-83.69988 


-83.45810 


-83.45810 


-83.45964 


{Eh) 


37.94086 


36.82427 


36.83193 


36.58714 


36.57918 


36.83133 


36.58747 


36.58747 






-4.86173 


-4.86641 


-4.86778 


-4.84360 


-4.84245 


-4.86771 


-4.84374 


-4.84373 




(Ec) 


-0.73575 


-0.73521 


-0.73530 


-0.73374 


-0.73366 


-0.73530 


-0.73374 


-0.73374 




(vr) 


32.34338 


31.22265 


31.22885 


31.00981 


31.00306 


31.22832 


31.01000 


31.01000 


31.01078 


{Enu) 


11.66667 


11.66667 


11.66667 


11.66667 


11.66667 


11.66667 


11.66667 


11.66667 


11.66667 


{V) 


-39.42376 


-40.79570 


-40.80503 


-40.78074 


-40.77317 


-40.80489 


-40.78144 


-40.78144 


-40.78219 


(Eel) 


-40.08293 


-41.28318 


-41.28437 


-41.37545 


-41.37535 


-41.28455 


-41.37566 


-41.37566 


-41.37566 


(E) 


-28.41626 


-29.61651 


-29.61770 


-29.70878 


-29.70868 


-29.61789 


-29.70900 


-29.70900 - 


-29.70899" 


N 


13.89834 


13.99939 


13.99865 


14.00002 


14.00003 


13.99864 


14.00000 


13.99999 


13.99998 


■"This 


is from the 


grid-DFT calculation; 


the corres 


ponding g 


rid-free DFT value is - 


-29.71530 a.u 





(ii) and (ill) 10 ^ a.u. seemed appropriate while for (i), we used both 10 ^ and 10 ^ a.u. 
(see Section III). 



III. RESULT AND DISCUSSION 

At first, we examine the convergence and stabihty of our nonrelativistic ground state 
total energy as well as all the individual energy components for CI2 molecule with respect 
to the grid spacing hr and number of grid points Nr {r & x,y, z) at a.n internuclear distance 
4.20 a.u. in Table I. We did a series of test calculations and here we present a select 8 of 
them. These are: (i) Sets A, B with hr = 0.3 and 0.4, both having A^,, = 32; (ii) Sets C, D, 
E with hr = 0.2,0.3,0.4, all having A^^ = 64; (iii) Sets F, G with K = 0.1,0.2, both having 
Nr = 128; and (iv) Set H with hr = 0.1, Nr = 256. All the calculations in this table are 
performed with an energy convergence criterion of 10"'' a.u. The reference values denote 
the results obtained from the GAMESS suite of quantum chemistry program j^^] with the 
same XC combination, basis set and effective core potential. The Hay-Wadt (HW) valence 



basis set 



25 



26| is used for all the calculations done in this work, where the valence orbital 



is essentially split into inner and outer components (described by two and one primitive 
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gaussian functions respectively). Expectation values for the following energy operators are 
reported: kinetic energy (T), total nucleus-electron potential energy (^"^), two-electron 
Coulomb repulsion energy {Eh), exchange energy {Ex), correlation energy {Ec), total two- 
electron potential energy (l^^*^), nuclear repulsion energy {Eny), total potential energy {V) 
((^"'^) + {^t^'^) + {Enu)), total electronic energy {E^i) and total energy {E) respectively. 
The Hartree, exchange and correlation energies of the reference values are not reported 
as these individual contributions were not available in the reference output. Furthermore, 
the total integrated electron density is given as A^, which is a good measure of the quality 
of the results. We note that the reference values [2^ for the total energy obtained from 
grid-DFT and grid-free DFT calculations are —29.70899 and —29.71530 a.u. respectively. 
The former uses the default "army grade" grid using Euler-MacLaurin quadratures for the 
radial integration and Gauss-Legendre quadrature for the angular integrations 11-3]. The 



grid-free calculations [l9|, l20| ideally use the resolution of identity to simplify the evaluation 



of molecular integrals over functionals, rather than the quadrature grids. The latter is quite 
appealing, for it enables one to avoid the finite grid and associated error; however this 
usually requires an extra auxiliary basis set to expand the identity. The first thing to note 
is that the maximum deviation in energy (off by almost 1.30 a.u.) from the reference value 
is shown by Set A; presumably the box length is insufficient to capture all the important 
contributions. This is also reflected in the poor value of A^. As the spacing is increased to 
0.4, the results for all the quantities get significantly better in Set B. As we further enlarge 
the box size with an increase in hr to 0.5 a.u. (not shown in the table due to lack of space), 
the total energy reaches a value of —29.71012 a.u., which is lower than the reference value 
by 0.00113 a.u. This is a reasonable agreement with the reference value; however there may 
be slight caused by the combined effect of a large h,. and small A^,.. It is noticed that Sets 
C and F both produce very similar results for all the quantities including A^, as one expects 
intuitively, for they cover the same box length. Set B also corresponds to the same box 
size as the above two sets and indeed the total energy is again comparable; however the 
individual energy components and A^ show slight deviations. In Set D, our results for all 
the quantities including A^ are very nicely matching with the reference values. Keeping A^^ 
fixed and increasing hr to 0.4 a. u. in Set E, the total energy deviates by 0.00010 a.u. from 
Set D and A^ remains almost unchanged. Thus among A-E, Sets D, E are to be considered 
two best values. In order to gain further confidence and examine the variationality, several 
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TABLE II: Variation of all the energy components and total number of electrons with respect to 
the grid parameters for HCl along with the reference values at R = 2.40 a.u. All quantities are in 



a.u. 





Nr = 32 




Nr = 64 




Nr 


= 128 


Ref. [271 


Set 


B 


C 


D 


E 


F 


G 




hr 


0.4 


0.2 


0.3 


0.4 


0.1 


0.2 




(T> 


D. 17590 


6.17910 


6.17900 


6.17510 


6.17909 


6.17796 


6.17811 


(vr) 


-37.16617 


-37.17230 


-37.17252 


-37.16470 


-37.17228 


-37.16998 


-37.17042 


(Eh) 


15.78969 


15.79289 


15.79266 


15.78825 


15.79287 


15.77980 






-2.74695 


-2.74742 


-2.74737 


-2.74674 


-2.74742 


-2.73626 






-0.41724 


-0.41727 


-0.41727 


-0.41723 


-0.41727 


-0.41706 




{vn 


12.62550 


12.62820 


12.62803 


12.62428 


12.62818 


12.62648 


12.62677 


(Enu) 


2.91667 


2.91667 


2.91667 


2.91667 


2.91667 


2.91667 


2.91667 


{V) 


-21.62401 


-21.62744 


-21.62554 


-21.62375 


-21.62743 


-21.62683 


-21.62699 


(Eel) 


-18.36478 


-18.36501 


-18.36551 


-18.36532 


-18.36501 


-18.36554 


-18.36555 


(E) 


-15.44811 


-15.44834 


-15.44884 


-15.44865 


-15.44834 


-15.44887 


-15.44889*^ 


N 


8.00023 


7.99992 


8.00002 


8.00003 


7.99992 


8.00000 


8.00003 



"This is from the grid-DFT calculation; the corresponding grid-free DFT value is —15.44888 a.u. 

extra calculations (Sets F, G and H) are done in a relatively larger and finer mesh. Thus 
Sets E, G, H all correspond to the same box length of 25.6 a.u. However, from E to G, total 
energy changes by only 0.00032 a.u., while G and H Set results are virtually identical. It 
is quite gratifying that Set G and H results not only completely match with each other in 
total energy; the component energies are identical up to 5th decimal place. It may also be 
mentioned that a calculation for A^^ = 128, = 0.3 (Set I) (not shown in the table) leads 
to an energy value of —29.70909 a.u., which again exactly coincides with Sets G and H. On 
the basis of above discussion, it is abundantly clear that Sets D, E, G, H are our four best 
results; while Sets D, E are sufficiently accurate for all practical purposes. 

Next in Table II we report results for the heteronuclear diatomic, HCl. Our basic pre- 
sentation strategy remains the same as in Table I, although in this case a fewer number 
of grid sets are given. Thus a total of 6 sets B-G are presented. As in Table I, here also 
we used an energy convergence criterion of 10~^ a.u. The reference total energies from the 
respective grid and grid-free calculations of —15.44889 and —15.44888 a.u. seem to be in a 
better agreement with each other than that for CI2. Overall, very similar trend is observed 
in this case as the grid parameters are changed, as for CI2. Note that Set B was inadequate 
for CI2; but for HCl this seems quite reasonable. All of these sets are capable of reproducing 



9 



TABLE III: Comparison of the calculated eigenvalues of CI2 and HCl with the reference values 
(grid-DFT). Negative values are given in a.u. 



MO 






Cl2{R= 


=4.2 a.u.) 






MO 






HC1{R= 


=2.4 a.u.) 






Set 


D 


E 


G 


I 


H 


Ref. [27] 




B 


C 


D 


F 


G 


Ref. [27] 


2<7g 


0.8266 


0.8274 


0.8268 


0.8268 


0.8268 


0.8267 


2ct 


0.7788 


0.7784 


0.7787 


0.7784 


0.7786 


0.7786 


2(T„ 


0.7167 


0.7175 


0.7168 


0.7169 


0.7168 


0.7168 


3(T 


0.4228 


0.4225 


0.4230 


0.4226 


0.4228 


0.4228 


So-g 


0.4332 


0.4340 


0.4334 


0.4334 


0.4334 


0.4333 


Iti 


0.2862 


0.2861 


0.2866 


0.2861 


0.2864 


0.2864 




0.3515 


0.3517 


0.3516 


0.3517 


0.3516 


0.3516 


iTTy 


0.2862 


0.2861 


0.2866 


0.2861 


0.2864 


0.2864 




0.3515 


0.3517 


0.3516 


0.3517 


0.3516 


0.3516 


















0.2859 


0.2861 


0.2861 


0.2861 


0.2861 


0.2860 


















0.2859 


0.2861 


0.2861 


0.2861 


0.2861 


0.2860 

















the reference values nicely up to the third place of decimal, although Set B result is, to 
some extent, inferior to the other five sets in terms of the individual energy components and 
N. All the energies, however, are above the reference value. Just as in previous table, a 
gradual improvement is observed as one passes from Set B-C-D. Moreover, Sets C and F 
results show mutual agreement with each other as in Table I. However for CI2, their total 
energies were above the reference value by roughly 0.091 a.u. and was correct only up to 
second decimal-place; for HCl, on the other hand, these are above the reference value by only 
0.00055 a.u. and A^ shows fourth place of decimal accuracy. As one passes from Set D to 
E, A^ remains almost unchanged and total energy increases by approximately 0.0002 a.u. (a 
trend also observed for CI2). Further calculations with A,, = 128, hr = 0.3 (Set I), produces 
exactly identical result as in Set G expectedly, and are thus not presented separately. Not 
surprisingly, as in CI2, in this case also our three best results are those from Sets D, E and 
G; with the former two being sufficiently accurate for all purposes, once again. 

Table HI lists the computed eigenvalues for several sets for CI2 and HCl at the respective 
R values as in Tables I and II, along with the reference values. For CI2, the calculated 
eigenvalues are either completely matching or show an absolute maximum deviation of only 
0.0001 a.u. for the Sets D, G, I and H; Set E gives an absolute maximum deviation of 
0.0007 a.u. (for 2ag, lo^ and Sa^). This is also apparent from a consideration of their 
performances in Table I. For HCl, Sets B, D, F give an absolute maximum deviation of 
0.0002 a.u., while the same for Set C is 0.0003 a.u. Set G results completely match with the 
reference eigenvalues. 

For further examination. Table IV displays the calculated negative total energies of CI2 



10 



TABLE IV: Comparison of the calculated potential energy curve of CI2 and HCl for several grid 



parameters along with the reference values (grid-DFT). Negative values are given in a.u. 



R (a.u.) 


CI2 (Total energy relative to — 2£ 


1 a.u.) R(a.u.) 


HCl(Total energy relative to 


— 15 a.u.) 


Set 


D 


E 


G 


I 


Ref. [27| 




B 


c 


D 


Ref. [27| 


3.50 


0.6508 


0.6504 


0.6509 


0.6508 


0.6509 


1.60 


0.1327 


0.1329 


0.1330 


0.1331 


3.60 


0.6d97 


0.6694 


0.6698 


u.DDy7 


n c t2c\o 

U.dd9o 


1.70 


r\ 00m 

0.2293 


0.2294 


U.2296 


0.2297 


3.70 


0.6ooy 


U.D806 


0.6840 


U-oooy 


0.6840 


1.80 


0.3004 


0.3005 


U.oUO/ 


0.3008 


3.80 


0.6943 


u.Dy4u 


U.D944 


0.6943 


u.by44 


1.90 


0.3520 


0.3523 


0.3525 


0.3526 


3.90 


0.7014 


0.7012 


0.7015 


0.7015 


0.7015 


2.00 


U.389U 


0.3893 


0.3895 


0.3896 


4.00 


0.7059 


0.7057 


0.7061 


U.7UDU 


0.7060 


2.10 


0.4147 


0.4150 


0.4152 


0.4153 


4.10 


0.7082 


0.7080 


0.7084 


0.7083 


0.7084 


2.20 


0.4317 


0.4321 


0.4324 


0.4325 


4.20 


0.7088 


0.7086 


0.7090 


0.7089 


0.7090 


2.30 


0.4420 


0.4425 


0.4430 


0.4432 


4.30 


0.7079 


0.7077 


0.7082 


0.7081 


0.7081 


2.40 


0.4481 


0.4483 


0.4488 


0.4489 


4.40 


0.7059 


0.7057 


0.7062 


0.7061 


0.7062 


2.50 


0.4501 


0.4505 


0.4508 


0.4509 


4.50 


0.7030 


0.7029 


0.7033 


0.7032 


0.7033 


2.60 


0.4494 


0.4498 


0.4501 


0.4502 


4.60 


0.6994 


0.6993 


0.6997 


0.6996 


0.6997 


2.70 


0.4466 


0.4472 


0.4473 


0.4474 


4.70 


0.6952 


0.6953 


0.6956 


0.6955 


0.6956 


2.80 


0A427 


0.4429 


0.4431 


0.4431 


4.80 


0.6906 


0.6908 


0.6911 


0.6910 


0.6911 


2.90 


0.4369 


0.4375 


0.4378 


0.4378 


4.90 


0.6858 


0.6861 


0.6864 


0.6863 


0.6864 


3.00 


0.4303 


0.4310 


0.4316 


0.4317 


5.00 


0.6809 


0.6812 


0.6816 


0.6815 


0.6816 


3.10 


0.4243 


0.4249 


0.4251 


0.4251 


(relative 


to -29 


a.u.) 


for Sets 


D, E, 


G and I {N., 


= 12J 


I, K = 


0.3) in 


columns 2- 


-5 with the 



reference values (column 6) covering a broad bond length region of 3.50-5.00 a.u. In the 
right panel, the same for HCl is given for three Sets B, C, D along with those obtained from 
the reference calculations (all relative to —15 a.u.) for R = 1.60 — 3.10 a.u., in columns 8-11. 
These are depicted in Fig. 1 for smaller R ranges to show the energy changes more clearly. 
We note that all these calculations are performed with a less stricter energy convergence 
criterion of 10~^ a.u. For CI2, all these four sets reproduce the qualitative shape of the 
potential energy curve very well for the entire range. Set D produces energy values quite 
well (higher by only 0.0001 a.u.) until R = 4.00 a.u., and thereafter develops a gradual 
tendency to deviate more. Nevertheless the maximum deviation is quite small (only 0.0007 
a.u.), that occurs for R = 5.00. The other two sets G and I are either completely matching 
with the reference values or show an absolute maximum deviation of 0.0001 a.u. We also note 
that the computed energy values have always remained above the reference values except 
for two instances [R = 4.00 and 4.30 for Set G). And leaving aside these two R values. Set 
G shows exact quantitative agreement with the reference curve. For HCl also, the three 
sets give very good qualitative agreement for the whole range of R as seen from their nearly 
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I 1 1 1 1 1 -0.448 I ' ' ' 1 

4 4.1 4.2 4.3 4.4 2.4 2.6 2.8 3 

R(a.u.) R(a.u.) 

FIG. 1: Potential energy curves for Cl2(left panel) and HCl(riglit panel) for different grid parame- 
ters. 

identical shapes. The absolute maximum discrepancies for the three sets B, C, D are 0.0014, 
0.0007 and 0.0002 for R = 3.00, 2.30 and 2.30 a.u. respectively. The best matching is 
observed with Set D. Obviously, as demonstrated for CI2 case, even more accurate results 
could be obtained with a fine-tuning of the grid parameters and refining the convergence 
criteria. 

Once the stability and reliability of our calculation is established, now Table V gives 
the computed kinetic, potential and total energies as well as for selected 10 molecules 
and 3 atoms (HCl and CI2 are omitted as they have been discussed earlier) to judge its 
applicability for a larger set of chemical systems. These are ordered in terms of increasing 
A^; corresponding grid-DFT values are quoted for comparison. All these calculations in 
this table are performed using grid Set E, which was found to be quite satisfactory for HCl 
and CI2. However, it may be mentioned that, for all these systems, using a smaller grid with 
Nj. = 32, hj. = 0.5, we obtained converged results of correspondingly similar accuracy as for 
CI2 (see discussion on Table I). For all of these species, we see there is excellent agreement 
with the reference values for all the quantities. In several occasions, energies are identical 
as the reference (for example As, Br2). The maximum deviation is observed for Na2Cl2 (by 
0.00056 a.u.). This once again demonstrates the faithfulness of the present results. 

Finally in Table VI, to gain further confidence, the calculated ionization energies, — enoMO 
(in a.u.) are compared with the grid-DFT result [27|; atomization energies are also compared 
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TABLE V: Kinetic ((T)), potential {{V)), total (E) energies and N for several molecules and 
atoms. All quantities in a.u. PW=Present Work. 



System (T> -{V) -{E) N 





r VV 


-tiei. |Z/j 


r VV 


rxei. Iz/J 


r VV 


r>„f [071 


r VV 


rxci. 12/] 


iNa2 


n 1 /I c^n'z 
U.i4oU/ 


u. i44yy 


U.O28U0 


U.OZ /yi 


U. 38292 


0.38292 


1.99999 


nnnnn 
2.U0U0U 


IN all 


U.ODyoi 


U.ODylz 


i.zy / iz 


i.zyoy ( 


n 'zo'zci 
U. ( z /oi 


U. / 2 /oo 


i.yyyyy 


nnnnc; 
2.UUUU0 


r 






8. / OOUl 


8. / o4U4 


O.ooO/O 


6.380 ( 1 


o.OUOUO 


/I nnnnn 

4.yyyyy 


As 


2.07461 


2.07354 


8.10154 


8.10047 


6.02693 


6.02693 


5.00000 


4.99999 


Br 


4.22038 


4.22011 


17.28157 


17.28131 


13.06119 


13.06120 


7.00000 


6.99967 


NaCl 


5.77569 


5.77639 


20.92133 


20.92235 


15.14564 


15.14596 


8.00001 


8.00059 


H2S 


4.90204 


4.90197 


16.10707 


16.10698 


11.20503 


11.20501 


8.00000 


7.99989 


PH3 


4.08953 


4.08953 


12.27387 


12.27383 


8.18434 


8.18430 


8.00000 


7.99965 


Br2 


8.55754 


8.55716 


34.74793 


34.74755 


29.19039 


29.19039 


14.00000 


14.00003 


H2S2 


8.75379 


8.75342 


30.00250 


30.00213 


21.24872 


21.24871 


13.99999 


13.99996 


MgCla 


11.62114 


11.62208 


42.34513 


42.34621 


30.72399 


30.72413 


16.00004 


15.99957 


Na2Cl2 


11.55066 


11.55242 


41.92870 


41.92990 


30.37804 


30.37748 


16.00002 


15.99686 


SiH2Cl2 


13.95036 


13.94989 


48.78729 


48.78685 


34.83693 


34.83696 


19.99999 


20.00015 
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34|. 
27|. 



with the experimental results [32| besides grid-DFT of [27| and other DFT results 
Our computed results for both these quantities are in excellent agreement with those of 
However, the atomization energies in several occasions show substantial discrepancy from 
experimental values and other DFT results. Note that the experimental values include zero- 
point vibrational corrections as well as relativistic effects. The largest error is found for CI2 
(about 13.5 kcals/mole). The results of 33l.l34| are based on all-electron calculations. Former 
used the 6-311++G(d,p) basis set and employed a combination of approximate exchange SC- 
a (recovering all the behavior of exact exchange) and a scaled GGA correlation functional. 
The latter used an optimum GG A/exact-exchange DFT. Probably use of more appropriate 
basis set and better XC functionals would further improve the present results. 

A few remarks may be made before passing. In this work our primary motivation was 
to demonstrate that the real-space CCG coupled with the Fourier convolution technique 
as employed here for the Hartree potential, could deliver accurate, physically meaningful 
results of "chemical" accuracy for molecular systems, through a small representative sets of 
12 molecules and 3 atoms. Thus no effort was made to reproduce either the most accurate 
theoretical results or the experimental values, which may be considered in future. These 
would inevitably require the inclusion of more extended and elaborate basis sets (containing 
the polarized and diffuse functions) as well as more accurate and sophisticated nonlocal XC 
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TABLE VI: Highest occupied molecular orbital energy, — eHOMo(a.u.) and atomization energies 
(kcals/mole) for molecules. PW=Present Work. 



System 




eHOMo(a.u.) 




Atomization 


energy (kcals / mol) 




PW 


Rcf. [27] 


PW 


Ref. [27| 


Other DFT 


Expt. [32] 


Na2 


0.1050 


0.1051 


14.64 


14.64 


16.15", 22.1' 


17.0 


NaH 


0.1458 


0.1457 


47.17 


47.17 




47.2 


HCl 


0.2863 


0.2864 


104.16 


104.16 


98.38", 101. 6*" 


102.2 


NaCl 


0.1811 


0.1810 


97.62 


97.63 


88 18" 96 3*" 


97.4 


HaS 


0.2267 


0.2265 


178.65 


178.65 


172.5' 


173.2 


PH3 


0.2323 


0.2323 


241.56 


241.55 


227.14", 227.1* 


227.1 


CI2 


0.2861 


0.2860 


43.72 


43.73 


40.18",59.0' 


57.2 


Bra 


0.2542 


0.2540 


42.68 


42.68 




45.4 


H2S2 


0.2371 


0.2370 


222.03 


222.04 




229.6 


MgClz 


0.2803 


0.2805 


185.14 


185.14 




187.4 


Na2Cl2 


0.2126 


0.2125 


248.33 


248.34 




243.1 


SiH2Cl2 


0.2903 


0.2905 


337.55 


337.55 




341.8 



^Ref. 
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34| 



functionals (with gradient and Laplacian corrections) having correct short- and long-range 
properties. At this stage, a few words on scahng is in order. Denoting A^";, and Ng as the 
number of grid points and basis functions respectively, one can assign the following scaling 
relations to the four most important steps of the whole process: (a) construction of the lo- 
calized basis set scales as Ng, (b) generation of the basis function in the grid scales as NbNg, 
(c) calculation of the electron density in the grid scales as N^Ng, and finally (d) build-up 
of the one- and two-body matrix elements of the Fock matrix scales as and N^Ng re- 
spectively. More detailed analysis of the scaling and computational time with respect to 
the system size may be considered in future communications. It is worthwhile mention- 
ing that our main motivation in this work is to build up a stable platform which enables 
us to perform the real-time dynamics studies (such as multi-photon ionization, high-order 
harmonic generation, etc.) of polyatomics in presence of an intense laser field that utilizes 
and exploits the enormous advances made in LCAO-GTO-based molecular DFT approaches 
over the years through many pioneering works. In taking up that, it was felt that TD im- 
plementation of the overwhelmingly successful ACG-based codes could be quite complicated 
and CCG-based approaches might be easier and straightforward to handle without making 
serious compromise on the accuracy and reliability. Thus the purpose is not to develop 
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another electronic structure code when there are several highly elegant and versatile quan- 
tum chemistry packages are available; but to make a modest base which provides sufficient 
accuracy and reliability to pursue the dynamical calculations, as mentioned. Although one 
could envisage some inaccuracies associated with the direct numerical integrations on the 
grid (due to the incompleteness of the grid), as our results suggest, with a proper adjustment 
of the grid parameters, this may not be so detrimental, at least for the systems of this work 
and use of non-uniform variable or adaptive grids may partly alleviate this problem. We 
note that we have also made some test calculations for systems hke P4, Sg (to be presented 
elsewhere), etc., leading to similar kind of results as the present work. This further vahdates 
the significance of this work. 



IV. FUTURE AND OUTLOOK 



We presented a detailed study on the performance of CCG in the context of atomic 
and molecular DFT calculations. The viabihty and feasibility of this is demonstrated by 
applying it to a set of 12 molecules and 3 atoms. This was achieved through an accurate 
representation of the classical Hartree potential through a Fourier convolution method on 
the real grid. Core electrons were represented by the HW pseudopotential and local LDA 
XC potentials were employed. An analysis on a cross-section of the calculated quantities 
including energy components, eigenvalues, ionization energies, atomization energies, as well 
as the potential energy curves, reveal that this can offer fairly accurate and reliable results. 
The results are variationally well-founded. More elaborate and systematic investigation 
on the properties such as atomization energies, vibrational properties, reaction energies, 
etc., of different chemical systems like small clusters, weakly-bonded molecules would be 
required to assess the success and performance of this approach. Better basis sets and 
XC functionals would be needed for further improvement. It may be interesting to study 
the performance of this in the context of all-electron calculations. Other interesting areas of 
study would constitute the real-time dynamics of small to medium size molecules in presence 
of a strong external TD field, or the chemical descriptor analysis through various local and 
global quantities (like softness, hardness, fukui function, etc). Some of these issues are 
currently being investigated by us. 
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